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Abstract 

We study the dependence on the external pressure P of the velocities vlt 
of long wavelength sound waves in a confined 2D h.c.p. lattice of 3D elastic 
frictional balls interacting via one-sided Hertz-Mindlin contact forces, whose 



T3 

a 

q ■ diameters exhibit mild dispersion. The presence of an underlying long range 

order enables us to build an effective medium description which incorporates 

>' 
CO ■ the radial fluctuations of the contact forces acting on a single site. Due to the 



non linearity of Hertz elasticity, self-consistency results in a highly non-linear 
differential equation for the "equation of state" linking the effective stiffness 
of the array with the applied pressure, from which sound velocities are then 
obtained. The results are in excellent agreement with existing experimental 

results and simulations in the high and intermediate pressure regimes. It 

O " 

O ■ emerges from the analysis that the departure of vl{P) from the ideal P 1 / 6 



Hertz behavior must be attributed primarily to the fluctuations of the stress 



field, rather than to the pressure dependence of the number of contacts. 
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I. INTRODUCTION 

Sound propagation in a dry confined granular medium still is, to a large extent, an open 
question. In particular, a long-standing puzzle is concerned with explaining the dependence 
of sound velocities on the externally applied pressure P. 

The load bearing intergrain contacts which ensure the mechanical stability of the packing 
are of the Hertz type, i.e. their longitudinal stiffness (along the intercenter axis) scales as 
P 1 ' 3 , with P the corresponding load. As shown by Mindlin |I|], this scaling also holds, for 
frictional balls, for their shear stiffness, provided that the shear load borne by the contact 
is much smaller than the friction threshold. So, one intuitively expects that the velocity of, 
say, longitudinal sound vl ~ P 1 / 6 . 

However, experimental results depart strongly from this expectation : vl is found to 
exhibit a much faster P-dependence, which is commonly characterized by "effective expo- 
nents" v = d(logf )/d(logP). Values of v of order roughly 1/4 are often mentioned. Such a 
behavior is observed for 3D random grain packings, which present what we will term " strong 
topological disorder" but also, more surprisingly, for artificially built regular arrays of quasi- 
monodisperse balls. This was first shown on a 3D f.c.c. lattice by Duffy and Mindlin 0, who 
found that the faster dependence with v ~ 1/4 in the intermediate pressure range tended 
asymptotically, at high P's, towards the Hertz P 1 / 6 dependence. Recently, Coste and Gilles 
|| have investigated in detail sound propagation in a 2D hexagonal lattice of steel balls. 
Their experimental results for vl{P), which are qualitatively similar to those of Duffy and 
Mindlin, have motivated the present study. 

Various explanations for this behavior for strongly disordered packings have been put 
forward. In particular, Goddard [|J proposed that a 1/4 value of v might originate from the 
existence of conical contacts, while de Gennes || suggested that the presence of heteroge- 
neous shells surrounding the grain bodies might in some cases be relevant. A possibly more 
natural explanation lies in the pressure dependence of the numberof load bearing contacts 
in the packing. It has been at the center of several recent works. In particular, a direct 



numerical study by Makse et al. |6j demonstrates the correlation between number of contacts 
and sound velocity in a 3D system. 

For regular arrays, such an explanation may at first appear doubtful. However, J.N. Roux 
has studied numerically the 2D h.c.p. structure and found that even a minute dispersion 
in ball radii leads to a similar effect. Namely, due to purely geometrical constraints, as 
the pressure increases, the average number of contacts per ball N c varies from ~ 2.5 at the 
rigidity threshold to its maximum value of 6 at high P. 

Although such effects seem to be within the reach of a mean field description, the various 
attempts in this direction have up to now failed to account for deviations from the Hertz 
power law ||, [|], [|K|. As already suggested by Makse et al. [|], we believe this to result 
from the central but implicit assumption of these models that the local contact arrangement 
deforms homothetically when P changes. This amounts to neglecting the essential effect of 
local stress inhomogeneities: due to the elastic Hertz deformations, they result in a dispersion 
of intercenter distances, and hence of the bond strengths. In other words, even in the absence 
of a change in contact number, the change of bond stiffness induced by a pressure change is 
non affine. 

While it is certainly very difficult to improve upon this approximation in the case of strong 

topological disorder, the case of a periodic array of weakly disperse balls seems amenable 

to a realistic description. Indeed, in this case, the existence of a reference lattice permits to 

formulate a mean field theory in the spirit of the single center self-consistent CPA approach 

IIJ , developed already long ago to describe the electronic properties of metallic alloys. 



Such a route was already explored in a series of articles on depleted elastic networks fL2| , 
13| . However, in these works, the distribution of active contacts was assumed to be known 



and independent of the external stress, and the self-consistent condition was formulated in 
terms of a single bond approximation, in the spirit of Kirkpatrick's approach to percolation 



problems fH}| . 



In this paper, we build an effective medium description of a 2D h.c.p. array of Hertz- 
Mindlin balls which does account for local deformation due to the disorder in ball radii. In 



contrast with previous theories, our self-consistency condition does depend on the global 
external stress. Clearly, in such a system, the higher the external pressure, the smaller the 
relative disorder. So, our mean field appears as a high-P expansion. It therefore comple- 
ments the numerical studies of Roux, which deal with the low-P regime where percolation 
effects are dominant. 

We show that our predictions account quantitatively for the experimental results of Gilles 
and Coste. Moreover, comparison with Roux's results in an overlapping intermediate pres- 
sure range allows us to determine the range of validity of our effective medium approach. We 
find that it holds down to pressures where N c has decreased by about 15% of its saturation 
value. From all this, we conclude that the basic physical effect responsible for departures 
from the P 1 / 6 law is, rather than the variation of contact number in itself, the disorder 
induced spatial stress fluctuations. 

The article is organized as follows. In Section II we set the basis of our model by writing 
the dynamical equations for a set of Hertz-Mindlin contacts under equilibrium forces aligned 
with the intercenter directions, and solve them for an ideal lattice of perfectly identical balls. 
In Section III, we build up our mean field description, apply it to the h.c.p. lattice, and 
obtain from it the equation of state, i.e. the force-displacement relation from which the 
effective bulk and shear moduli are derived. Section IV compares in detail the mean field 
predictions with the experimental and numerical results. 

II. BASIC MODEL AND DYNAMICAL EQUATIONS 

A. Equations of motion 

Let us first consider two spherical balls, labeled (i), (j), of radii di, dj, made of the same 
material of Young modulus E, Poisson ratio a, density p. The balls at equilibrium are in 
Hertz contact [|IJ under a force / directed along the intercenter axis ij, i.e. normal to the 
contact circle. 



The normal force / displaces the intercenter distance from dij = \(di + dj) to Oy = 
dij — dij, and: 

Due to solid friction, when submitted to a tangential force smaller than the friction 
threshold fi s f, the contact is pinned and cannot slide. The elastic response of this system to 
small additional forces Sf x , Sf y in the (x, y) plane is described, in the linear approximation, 
by two stiffness coefficients, given by the Hertz-Mindlin expressions, namely 0: 

(i) compression (normal) stiffness 

h 3 = -^ = (^* 2 ^/) 1/3 (2) 

where 



E >»\-\ _ 1 r J -i /-i. 



^=2(13^2) W- ^ + d ") ( 3 ) 

(ii) shear (tangential) stiffness : 

«<i = qki ^ = 2 _ ° o (4) 

Let us insist here that the fact that the ratio rj = Kij/kij is a mere material parameter, 
independent of the values of the equilibrium forces, only holds for the case of a purely normal 
equilibrium loading, which we assume here to be the case. In the general situation with a 
finite equilibrium tangential load f t , i] —* r][l — (ft/l^sf) 1 ^ 3 }- I n & U the following, we will 
restrict ourselves to situations where all (equilibrium and non-equilibrium) interball forces 
are near normal, so that both the static and dynamic shear responses are described by the 
linear shear stiffness (equation (|4])). 

Small displacements of the balls can be decomposed into rigid translations u i; Uj, and 
rotations about Oz by angles <fii,4>j. Denoting ny,ty the unit vectors normal and tangent 



to the (ij) contact with n^ directed from (i) to (J) (see Figure |l|), one immediately finds 
that the force and torque on ball (i) associated with the (ij) contact read respectively: 

SFij = kijKiii - \ij) .h^hij 

+7]kij[(u i -Uj).ii j ]ii j (5) 

-r}kij\{dj^)j + di(f>i)tij 

SQj = r]ki~[(uj - Ui).t y - -{dj<j>j + dkfa)] (6) 

and, in the small displacement limit appropriate to sound propagation, the equations of 
motion for a 2D lattice of balls are: 

M 4 u 4 = 5>F M M l = ^f- (7) 

J^ = E^m h = ^ (8) 

{j} 

where the sums are restricted to nearest neighbors in direct contact with (i). 

B. Ideal h.c.p. lattice 

We now consider the ideal case of a 2D h.c.p. lattice of balls of equal diameters d prepared 
so that, at equilibrium, the interball forces, of magnitude /, are directed along the normals 
to the contacts fij (see Figure [I]). Such a "hydrostatic" configuration can be realized by 
applying a force per unit length P = fy/3/d on a hexagonal container with walls along the 
dense ball rows, as realized in ||. The unit cell is defined by the two vectors ai j2 = dhi^, 
the ball centers by R mn = mai + na 2 . 

One then obtains from equations (H) - @, for the vibrational modes of this system: 

u mn = ue* (q - R — ^ <p mn = 0e l(q - R --^ (9) 

Mu 2 u = 

3 A „ 

2k E (1 - C P (q) [(u.n p )n p + r]u.t p )tp] (10) 

p=i 

3 

+ir)kd(j)J2 S P (c[)t p 

P =i 



Iu 2 <f>=-ir]kdJ2S p (q)(u.i p ) + ^-«^(1 + C p (q)) (11) 

p=i A p=i 

where k is the normal stiffness common to all contacts, and: 

S p (q) = sm[d(q.n p )] C p (q) = cos[d(q.n p )] (12) 

6 

and use has been made of relations such as J2 % = 0. 

P =i 

The exact spectrum in the full Brillouin zone, computed for q along two directions of 
high symmetry, is shown on Figure ^|. It can be inferred from these results that, in this 
close-packed lattice, the anisotropy of the spectrum is small. 

In the elastic continuum, long wavelength limit qd <C 1, decomposing the translation 
amplitude into its longitudinal and transverse components : u = Uj|q + u±(z A q), one finds 
that the vibration spectrum, which is isotropic due to the hexagonal symmetry of the lattice, 
is composed of three branches: 

(i) pure longitudinal acoustic modes of frequency ujt = viq, where the longitudinal sound 
velocity vl is given by: 

w|=-(l + ^)-^ F d 2 (13) 

L 8 V 3 J M v ; 

(ii) Two branches of mixed modes containing both a transverse translational and a rota- 
tional component. One of them is acoustic : u = Vxq, where the transverse sound 
velocity reads: 

4 = ^(1 + ^)^ 2 (14) 

The second, which corresponds to pure rotation in the q = limit, is an optical branch, 
defined by: 

^ = 30^(1-1^) (15) 



That is, as already shown by Schwartz et al. |T5| , the specificity of the vibration spectra 
of granular systems in frictional Hertz contact, as compared with atomic systems, lies in 
the additional degree of freedom associated with ball rotation. This remains coupled, in the 
long wavelength limit, with shear deformation, leading to a contribution (—3r]kd 2 /AM) to 
v\. 

Also shown on Figure ||] are the long wavelength dispersion curves. They can be seeen to 
provide a good approximation for the exact spectrum in a sizeable fraction of the Brillouin 
zone. 

It is worth recalling, at this point, that our equations of motion cease to be valid when the 
frequency approaches that of the lowest acoustic resonance of a ball : o; RES ~ v BVhK /d, with 
v bulk a sound velocity of the material constituting the balls. Indeed, in this situation, elastic 
deformations are no longer localized in a small region of extension on the order of the contact 
radius, internal deformations become important, and the restoring forces can no longer be 
described simply via the Hertz-Mindlin stiffnesses. Roughly speaking, this means that our 
expressions for the acoustic branches of the spectrum are valid provided that vl,t "C ^bulk- 
In view of equations flT3|), (|TJ|) and (g), this simply amounts to (f /Ed 2 ) 1 / 6 <C 1, which is 



realized under ordinary experimental conditions [|16| . Note that this condition is equivalent 
to stating that the radius of the Hertz contact circle must be much smaller that the ball 
radius, which is precisely the condition for the Hertz approach to hold. 

From planar continuum elasticity applied to a medium with hexagonal symmetry, the 
velocities of sound waves with propagation and polarization directions in the basal plane 
read : 

2 K+G 2 G 

Vt = Vtn = — (16) 

p p 

with p the mass density of the medium, K and G its bulk and shear moduli. 

Comparison between these expressions and equations fll^D , ( fUD enables us to define 
elastic moduli for our ball lattice - a result which will be of use in the disordered case. They 
read : 

8 



where p is related to the density of the ball material by: p = np/3\/3. 



III. DISORDERED LATTICE 

A. Random h.c.p. array of balls under hydrostatic compression 

Consider an ensemble of balls whose material parameters are identical, while their di- 
ameters vary at random with a continuous or discrete statistical distribution. A diameter 
value d®, with a formal label Q, has the probability, or fractional concentration, c^, and 

Ec° = i. 

The mean diameter d of a set of iV ^> 1 such balls may be obtained by configuration 
average, which we will denote by (...): 

d = T7 E di — > d = (d) = J2 c Q d Q (18) 

iV i Q 

We define the random deviations by 

Ai = di - d, A® = d Q -d 

(19) 

(A) = 0, Al us = (A*) = Ec°(^) 2 

In actual computations, we will use the uniform distribution with full width W = 
m&x{A®} — min{Z\^}, which means A RMS = W\f\2. 

One sample (configuration) of the random h.c.p. array of balls may be created by ran- 
domly distributing the balls in an uncorrelated way over the sites of a reference h.c.p. lattice. 
The number of balls is assumed to be sufficiently large for the thermodynamic limit to be 
approached. The sample is then random, but macroscopically homogeneous. The balls 
are brought into contact and further compressed by external compressive forces applied to 
the sample boundaries so that the average internal stress is hydrostatic. An easy way of 
achieving this is to assume a hexagon shaped sample and to apply to all its sides the same 
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macroscopic pressure force. The pressure force per unit length P will be denoted henceforth 
linear pressure . Under this pressure, the size of the compressed sample is reduced, while 
the symmetry of the lattice is preserved on average. It is thus meaningful to introduce an 
average lattice spacing a as a macroscopic parameter having a thermodynamic limit and 
globally characterizing the state of the sample. 

The pressure and the size of the compressed system are related by the macroscopic 
equation of state for the random lattice under hydrostatic compression. It will be convenient 
to introduce it in the form 



f = f(a), / = P-d/V(3) (20) 

where / is the functional dependence in question and / denotes the average intergrain 
hydrostatic force associated with the linear pressure P. The effective normal stiffness k, and 
the effective bulk modulus K EFF are then given by 



d/W K , o^_V3 - 

da ' ****'*- fda*\ ' 2 

dP 



fc (a) = _^Z ? Km . d= -- r ^ v = Z?..~ k (21) 



As for the d factor on the left hand side, cf. Eq. (|Hj). 

While the global characteristics, P and a, have thermodynamic limits, the ball positions 
and contact forces are subject to pronounced fluctuations, in particular for small external 
loads. In experiments, this double nature of the disordered state is manifested by the 
coexistence of coherent signals and irregular speckles in the acoustic response 0, fl7f 
numerical simulations most clearly by the formation of force chains [[Fl 



in 



Let us first verify that the geometrical disorder can be taken as small, as demanded is 
Sec. IIA. If the relative diameter spread is small, the equilibrium disordered system can be 
treated, for any external pressure, as a distorted lattice. The balls are slightly displaced, 
and, in the case of frictional balls, the contact points may be somewhat off the intercenter 
lines, which gives rise to non- normal contact forces. Taking as representative the data of [|3] , 
W < 8/im, d = 8 mm for steel balls, we see that the relative dispersion in "bond" lengths 
is of the order of 10 -3 . By geometrical considerations, this corresponds to deviations of 

10 



the direction of the contact force from the normal also of about 10~ 3 radian. These figures 
indicate that the geometrical disorder in this case is small indeed. 

Fluctuations in the magnitude of the random contact forces, by contrast, may be quite 
large. By Eq. ([!]), the contact force as a function of intercenter distance is given by 



fij{x) 



§£*(d£)s(dy-X 



d 







d{j — x > 



dij — x < 



(22) 



As indicated by the underbrace in (P^), it is consistent to neglect the randomness of d*-. 
Namely, the fractional fluctuation of d*- ~ d + \Ai + Aj is ~ A RMS /d. 

On the other hand, the disorder of the last factor in fl2"2| ) is crucial: for the Hertz 
displacement dy — x, the ratio A RMS /(d — x) may be large or small depending on the degree 
of compression of the balls, while both basic conditions, A RMS <C d (small disorder), and 
d — x <C d (Hertzian picture), remain satisfied. We may thus envisage three different 
regimes: 



regime 


condi tion 


low pressure 


^RMS -^ ™> •£ 


intermediate pressure 


^RMS ~ U X 


high pressure 


A nm < d-x 



In the high pressure Hertzian regime, A RMS < rf - i C d, the disorder appears as a 
perturbation of a basic homogeneous and homogeneously compressed crystal. This natural 
conjecture has given rise to two simple, but useful approximations (see [0]) for the true 
effective medium: 



(i) The Averaged Lattice Approximation (ALA) replaces the true sample by an ideal lattice 
of balls with the diameter d. The constitutive law and the bulk modulus thus read 



[0 



f(a) = \E*d*{d 



a) 2 



F(a) 



(23) 
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/ 9 P \3 

K ALA = E*(— —-) (24) 

Thus, the bulk modulus naturally obeys the plain Hertz | law, as is appropriate for a 
periodic array (see Sec. II.). 

(ii) The Averaged Force Approximation (AFA) assumes the balls to be positioned at any 
pressure exactly at the lattice sites, while the contact forces are given by (|22|). In this 
strictly symmetric geometry, no shear forces between the grains occur, and the effective 
hydrostatic force between two adjacent sites is obtained by configuration averaging, 

/(«) = (fij( a )) = F AVE (a) 

(25) 
F AVE (a) = £ c^c^F^'ia) 
Q,Q' 

where we have introduced the notation: 

FW(x) = lE*d^(d QQ ' - x)l (25a) 



Asymptotically, 



for a force at a separation x between two balls of the prescribed species Q, Q' with 



^~ F ^{ 1 + l(^hf) < 26 > 



and the bulk modulus becomes 



2 . f*JP*2j-l\i 



1 A* ■ (%E**d 



K AFA = Kjoj, ■ [ 1 - i- RMS K3 pl - ' 1 (27) 

These two equations show that the AFA effective contacts are non-Hertzian, and that 
the pressure dependence of the elastic modulus deviates from the simple | power law. 

While the ALA has an arbitrary nature, the AFA seems to be justified for high pressures, 
when the geometrical disorder is small compared to the compressive displacements, while 
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the contact forces still continue to exhibit random fluctuations. These force fluctuations 



then appear as being responsible for the non-hertzian features of the effective contacts |20 

At lower, "intermediate", pressures, according to our above classification, all grains are 
still mutually engaged through their contacts, but the random displacements become com- 
parable to the global compressive deformations. The AFA is then not sufficient, while an 
approximation taking these local lattice distortions into account in an averaged manner may 
be satisfactory. The EMA described in the next section is an attempt in this direction. 

B. Effective Medium Approximation 



Now, we are ready to build our Effective Medium Approximation for the disordered 
hexagonal array of balls. We use this name, conventional in the context of granular assem- 
blies, but, as sketched already in the Introduction, a more descriptive name would refer to 
the mean- field nature of the approximation, or, alternatively, to its " single site" character. 
The universal idea of such approximations is as follows JTTJ] . We will only consider disordered 
systems in which a periodic geometrical lattice (2D h.c.p. in our case) is filled by elementary 
objects associated with individual sites (balls for us) and having randomly variable char- 
acteristics (radii). It is assumed that the macroscopic properties of such random systems 
can be obtained, in principle, by configuration averaging. The configurationally averaged 
system is exactly periodic again. In the mean-field approximation, one assumes that it can 
be represented as a periodic array of effective elementary objects similar in nature to the 
random elements of the real array. The characteristics of these are determined by the follow- 
ing self-consistent procedure. One of the effective elements is replaced back ("substituted") 
by a true random one. The new system is locally sampled in a random fashion. It is then 
required that a configuration average of these locally distorted systems restores the average 
behavior. This condition determines the characteristics of the effective elements. Once this 
has been done, the whole task of configuration averaging is complete. 

An effective medium theory along these lines was developed by Feng et al. in |T2j and 
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|T3| to study depleted elastic networks on lattices. The basic element in the theory of these 
authors was a single bond connecting two lattice sites. The bonds were described by their 
linear stiffnesses, whose random distribution was prescribed a priori. The resulting EMA (a 
"single-bond Coherent Potential Approximation") provided a theory for the effective linear 
elasticity of the network, and for its vibrational spectrum. 

While this bond- EMA was a successful theory in its own area, we contend that it cannot 
be used for random Hertz lattices for two essential reasons: 

(i) The elementary object in a granular system is a ball. The star of (six for 2D 
h.c.p.) contacts surrounding the ball, that is of "bonds" stemming from the center, is 
statistically correlated, and the bonds cannot be treated as independent. 

(ii) The linear stiffnesses of the contacts are not known beforehand. They are indirectly 
specified by the average external pressure, but their local fluctuations depend on the 
equilibrium ball positions and cannot be determined independently of the non-linear 
static equilibration of the Hertz array at a given pressure. 

These two points are of a different nature and importance. The second point holds for any 
granular system, and we believe that it is precisely this which has been the obstacle against 
developing a satisfactory EMA for the acoustic response of granular materials. It may be 
said that the grain network should fulfill two contradicting roles at the same time, namely it 
should constitute the medium for wave propagation, and, as well, the random scattering field. 
We will see how this basic problem can be overcome in our rather specifically constructed 
case. 

1. EMA for frictionless balls 

To develop an EMA incorporating these features, we consider first the case of non- 
frictional Hertz contacts (formally, r\ — 0). The averaged array is assumed to consist of 
effective balls, whose diameter is a. The principal assumption is that the average hydrostatic 

14 



force f(a) can be interpreted locally as a contact force between two effective balls, so that 
k(a) = — 4-f{a>) is the stiffness of an effective contact between two such balls. Now, we 
select one site, "0", as central and substitute a ball with diameter dP for the effective ball. 
The difference dfl — a in diameters will give rise to an elastic deformation of the effective 
lattice. In other words, the substituted ball acts as an elastic inclusion. On average, the 
deformation should cancel. This would define the EMA condition, if only we were able 
to determine the force between a real and an effective ball. This is not possible directly, 
and we propose to overcome this by considering a cluster consisting of the central ball and 
a "corona" of its nearest neighbors, as sketched in Fig.||]. These neighbors have a hybrid 
nature. In other words, the interface between the effective medium and the inclusion passes 
through the corona balls. All contacts between the corona balls and the adjacent effective 
balls are taken as effective ones, associated with the force /. This type of approximation is 
based on the mean field reasoning: the fluctuations of the remote parts of the lattice should 
not have a significant effect on the central site. 

Inside the cluster, the corona balls might appear as "true" randomly chosen balls, so 
that the contact forces would be of the form F®® as given by Eq.(25a). A straightforward 
procedure involving averaging over the individually equilibrated positions of such randomly 
composed clusters of 7 balls, while conceivable, would be disproportionately clumsy. 

We prefer to introduce a model of the corona, which, while simple and transparent, 
captures the core features of the problem: 

• The contact forces between the central ball and the corona balls are averaged over the 
corona configurations: 

foi(x) -> F Q (x) = J2c Q 'F QQ '(x), 1 = 1-^-6 (28) 

Q' 

These forces incorporate in full the symmetric, radial fluctuations caused by the ran- 
domness of the central ball, while all angular correlations are averaged over and the 
forces have the full hexagonal symmetry. 
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• The contact forces between the touching corona balls are assumed to be given by the 
effective force law f(x). 

The displacements of the corona balls are then also symmetric. This restores the basic 
picture of a symmetric inclusion in the effective lattice. The central ball remains at its 
site, while the corona breathes around it symmetrically in accordance with the central site 
occupancy, and transfers this to an equilibrium symmetric distortion of the surrounding 
effective lattice. 

Returning to Fig. [| we see that the cluster is surrounded by 12 effective balls forming a 
hexagon. They are of two kinds : six occupy the corners ("on top" positions with respect to 
the corona balls), and another six sit along the edges ("bridge" positions). The displacement 
field around the "1" ligand is sketched in Fig. |j. The scalar u is the radial displacement 
common to all corona balls, u a and up correspond to the on-top and bridge neighbors, 
respectively. The arrows indicate the directions of the displacements, whose magnitude 
varies with the Q label of the central ball, and should be labelled vP, u x , . . ., etc.. 

Ball 1 is acted upon by contact forces from all its neigbors. The force coming from the 
central ball is random (Eq. ( p8| ) ) : 



ioi 



F Q (a + u Q ) fii (29) 



In the regime of small displacements (~ high pressures), we may expand 

foi = (F Q (a) - K Q (a)u Q ) n 1; K Q (x) = --^-F Q (x) (30) 

The remaining five forces are given by the effective interaction /. They are of three types. 
In the regime of small displacements, we obtain (cf. Eq. (^)): 

fai = /(a)n 4 - k(a)(u Q - m^) n 4 
f/3i = /(a)n 5 -k(a)- 

{{v9 n x - u Q p±{h x + n 2 )} • n 2 ) n 2 (31) 

f2i = f{o) n 3 - k(a)u Q n 3 
etc. 
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The forces depend on the three displacements u®,u®,Up. In the linear elasticity regime 
corresponding to our assumption of small displacements, a simple universal relation exists 
between these three displacements: 



u, 



au, up = (3u 



(32) 



where a and (3 are geometrical parameters independent of the magnitude of u. They can 
be determined for a linear elastic array of non-frictional balls once for ever, although not in 
a closed analytic form. Details of the calculation of both parameters are in the Appendix. 
The resulting values are 



a 


0.585405 


(3 


0.232971 



Now, we require that the total equilibrium force on ball 1 vanishes, and subtract from this 
condition the equilibrium condition for the effective lattice: 



(33) 



f i + f 6 l + f/31 + f«l + f/3'1 + f 2 l = 

/(a) • {iii + n 2 + n 3 + n 4 + n 5 + n 6 } = 

Only the component parallel to the 01 connecting line is non-trivial. With the help of Eq. 
(pl|) it reduces to a single equation for a single unknown, the corona radial displacement, 
with the solution 

F Q (a)-f(a) 



u 



Q 



KQ(a) + Ek(a) 



(34) 



where 



2 a P 2 



(35) 



is again a geometrical factor. The equation has a simple interpretation. For a given Q, 
the difference between the force exerted by the inclusion and the average force gives rise to 
a radial displacement. Its magnitude is related to the force by an additively renormalized 
stiffness reflecting the fact that the corona is supported by the rest of the (effective) lattice. 
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In the mean field context, this renormalization is a local field correction, and the parameter 
H measures its dimensionless strength. 

Now, we may impose the self- consistency condition, requiring that the average displace- 
ment be zero: 



w 



F Q -f 



>=5 c Wm=° (36) 

where the argument of all functions is a. This equation can be rearranged by introducing: 

^ave = (FQ)=Ec Q FQ 
1 Q l (37) 

The first quantity is easily recognized as identical with the averaged force of Eq. fl25|) of the 
AFA in the preceding subsection. The second quantity is an averaged radial stiffness of the 
lattice surrounding a single ball inclusion. The EMA condition (|35|) is equivalent to 

The average hydrostatic force for a given lattice spacing is seen to contain two contributions. 
First, the averaged contact force, and, second, a term having the characteristic structure of a 
second-order correlator, known from CPA-like theories. The second term is also responsible 
for the peculiar self-consistent character of the equation (|38|). It is not the unknown / itself 
which appears on its r.h.s., but its derivative k(a) = ~-^f( a ), an d (§^) has the form 

/(a) = F AVE (a) + F{Z \ a, j-J(a)) (39) 

where T is a complicated but well defined function. The EMA equation is thus a first 
order differential equation determining / as a function of a. In other words, we meet here 
a case where self-consistency cannot be written directly for an isolated value of the control 
parameter a, but involves intrinsically the whole functional dependence f(a). 

To select the physically relevant solution of (|3"9"D, we choose the high pressure asymptotic 
boundary condition 



f(a) ~ F AVE (a) for Z\ RMS < d - a ( < d) (40) 



It can be explicitly verified from(p8|) that the contact of / with F AVE stipulated by (|40f) is of 
second order, so that the boundary condition is formally justified. Physically, this condition 
means that the EM A coincides with the AFA in the asymptotic limit of high pressures, but 
extends the mean field description of the system to an intermediate pressure region. We may 
also say that our EMA is a self-consistent resummation of the high pressure perturbation 
expansion for the equation of state f{a). 

The solution of the EMA equation ( p9|) is obtained without difficulty by numerical iter- 
ation starting from f^(a) = F AVE (a). 

2. EMA for frictional balls 

The case of frictional balls is very important for a proper comparison with existing 
experimental data. A formal development of the EMA for frictional balls can proceed 
quickly, because the steps are almost identical with those made for the non-frictional case. 
It is necessary, however, to extend the description of the effective medium. As before, it 
is assumed to be composed of effective balls with diameter a, interacting through effective 
contacts. In addition to the average hydrostatic force / and the normal stiffness k = 
—4-f{a), we have to postulate also a local shear (tangential) stiffness R as a third basic 
characteristic of the effective contact. For R, we formulate our basic conjecture: 

R = r]k (41) 

In words, the effective shear and compressive stiffnesses are related in the same manner 
as the two stiffnesses for actual Hertz-Mindlin contacts (see Eq.(^)). We will discuss this 
intuitively appealing conjecture below, but first derive the EMA equations following the 
preceding paragraph step by step. 

The central ball inclusion and its symmetrized corona are introduced without changes, 
and both Figures ^| and ^ remain valid. Because of the symmetry of the cluster, there is no 
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torque acting on either the central or the corona balls, and only the force equilibrium has 
to be considered. The central ball is equilibrated automatically, and we inspect the forces 
acting on the corona ball 1, linearized again with respect to the small displacements. The 
full expression (|5|) must be used now. The relative displacements remain radial for all the 
neighbors of ball 1, except for the j3 and j3' ones (see Figure f|. For /3s, the contact force has 
a shear component, 

f/3i = /(a)n 5 

-~k(a) ({«« h, - u^(fti + n 2 )} • Aa) A 2 (42) 

-k(a) ({u Q ni - «^(ni + n 2 )} ■ t 2 ) t 2 

and similarly for (3' . Finally, the displacements u®,uQ,ui obey, in the linear elasticity 
regime, universal relations similar to 



u a = a(rj) u, up = (5{rf) u (43) 

where the a and /3 parameters depend now on the stiffness ratio r], as indicated. 

Introducing Eqs. ( 41|) - ( ^3|) into the equilibrium condition (0), which is fully general, 
we obtain the corona radial displacement as 

uQ = gw - /(„) 

KQ(a)+E v k{a) 



in complete analogy with the frictionless case, Eq. fl34|). The only change is an ^-dependent 
geometrical factor, which reads now 

S, = |-a(r?)-^)f 

(45) 

+ ^(|-/5(^)f) 
The EMA condition (tfi) = 0, Eq. (^), has a universal character, and with tft given 
by (pH]), it can be brought successively to various explicit forms following strictly Eqs.(p6|) 
- (|39|), with H of Eq.(^5|) replaced everywhere by E v as given by (|45|) . For convenience, we 
present here the final form of the EMA, 

/(a) = F AVE (a) + ^(S, | a, £/(a)) (46) 
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The corresponding boundary condition ( f40|) is ^-independent. The two complementary cases, 
frictional and non-frictional, are thus given by the same differential equation with the same 
boundary condition. They are only distinguished by the value of the E v parameter, which 
corresponds to 77 ^ 0, and 77 = 0, respectively. In other words, the difference between the two 
physical situations is reflected in the EMA through the strength of the local field correction 
appearing in the renormalized contact stiffness. 

C. Macroscopic properties in the EMA 

To complete the EMA analysis of our granular system, let us now derive the expressions 
for the sound velocities. Macroscopically, the averaged system is a periodic h.c.p. structure, 
so that it is acoustically isotropic and there exist just two sound velocities vl and vt, as in 
the ideal lattice of Sec. IIB. There, we proceeded from the full dispersion law in the Brillouin 
zone to the long wave limit (Eqs.([T3|), flli])), then to the elastic moduli (|T7|) on the basis of 
the macroscopic (continuum) relations flip]). 

For the disordered system, we shall go in the reverse direction. The EMA yields directly 
the equation of state f(a) and the effective compressive stiffness k(a). The equation of state 
is a monotonous function of a and thus can be inverted into a(f). Using Eq. (pD|)), we 
obtain finally the average lattice spacing as a function of the linear pressure P, and any 
function of a, such as the bulk modulus, can be converted into a function of the pressure. 



The bulk modulus is obtained from the stiffness k using Eq. ([H]). It should be stressed 
that Eqs. (|T7|) and ([H]) have an identical form; in fact, (|17D as far as it concerns the bulk 
modulus is but a special case of (|2l|) for zero disorder. As concerns the shear modulus, since 
our method only allows for treating the effect of hydrostatic stresses, we resort to an ansatz, 
i.e., we propose to extend relation ([17]) for G to the disordered case. We then get : 

k _^~ k r - Q + g} K IATJ\ 

-f^EFF — „ ~7; '-Jeff — _ -"-eff \*' ) 

While the first relation is exact, we should discuss the meaning of the second one. It 
turns out that the equivalent relation ([Ij]) for the ideal array is a direct consequence of the 
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Mindlin ratio (||) between near-normal stiffnesses and does not depend on the specific value 
or character of the compressive stiffness. Thus, the second of Eqs. (|£7|) is in fact equivalent 
to our conjecture ([H|) about the ratio of effective stiffnesses, k = rj k. Once the effective 
moduli are known, the sound velocities become 

2 -^EFF ~T (j"EFF 9 ^EFF ;, n i 

v L = ; , v T = — ^- (48) 



Equations ( (47)) then predict for the velocity ratio that 

v\ 3 + r] 



(49) 



v \ I + 77 

This is a very strong prediction, according to which the velocity ratio should depend neither 
on the pressure nor on the disorder level, being fully determined by the material properties 
of the grains. While conjecture ([41]) is not susceptible of direct verification, the predicted 
sound velocity ratio can be tested for a wide range of systems. Note that a similar conjecture 



concerning this ratio was formulated earlier by Schwartz and Johnson [ITH] 



IV. DISCUSSION 

We first confront our mean field predictions with the experimental data obtained by 
Gilles and Coste ||, obtained on an h.c.p. lattice of steel balls confined by a hexagonal 
frame under hydrostatic loading. They measured the time of flight of a low frequency sound 
pulse between two opposite sides of the frame. From this we obtain the longitudinal sound 
velocity t>L(/), where / is the average equilibrium contact force (Eq . (|20|) . 

The experimental data show that, at the highest loadings used, vl approaches closely 
the P 1 / 6 Hertz behavior. This confirms that, in this regime, the applied pressure is high 
enough for disorder effects to be weak, and for the compressed lattice to be close to the ideal 
h.c.p. one. 

We therefore focus first on the data for the largest forces used in the experiments. For 
example, an external force of ~ 990 N applied on each side of a hexagon with 31 balls in 
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contact with each side corresponds to / = 18.446 N. Taking for steel E = 9.20 10 10 Pa, 
a = 0.276, p = 7840 kg.m -3 , ball diameters d = 8mm, we get from Eqs. (0), (^) and 
(O) vl = 791m. sec -1 , to be compared with the experimental value of 778 m. sec -1 . The 



agreement is seen to be excellent, to the precision of the numerical input. 

Note that the contribution of the tangential contact stiffness is non negligible : with the 
above value of er, the Mindlin coefficient r] = 0.84. Neglecting the frictional character of 
the contacts (77 = 0) leads to underestimating the longitudinal velocity by ~ 11%, well out 
of the experimental error bar. For the transverse velocity, the two values would differ by 
~ 35%, but no experimental data are available. These differences in velocities are clearly 
reflected in the different slopes of the dispersion curves for the frictional and non-frictional 
cases in Fig. ||. 

With the values of vl at hand, we may return to the discussion at the end of Sec. 
IIB and check first the validity of the Hertz approach. The bulk sound velocities for the 
above material parameters of steel are t>L iB uLK = 5825m.sec -1 , wt.bulk = 3252m.sec -1 and 
we may conclude that the condition vl,t ^ v l,t,bvlk is obeyed. To verify the validity of 
the continuum limit, we have to invoke the ground frequency of the pulses used in the 
experiments, uo = 2tt x 6 500 Hz. This yields q ■ (2ir/d)~ l p» 0.066, much less than 1 indeed. 
Expressed in terms of wavelengths, A ~ 0.122 m ~ 15d. These estimates made for one 
pressure may be taken as representative for the whole experimental pressure range, as the 
experimental velocities vary between 500 and 800 m. sec -1 . 

In order to analyze the pressure dependence, we deduce from the vl data the effective 
normal stiffness k with the help of equations fl4*9|), (fi8|), assuming frictional contacts with 
7] = 0.84. It is plotted on Figure [5[ It is clearly seen that the high pressure data fit 
the straight line corresponding to the ideal lattice. As / decreases, the logarithmic slope 
increases, coming close to the popular v ~ 1/4 value. However, no sharp crossover is 
identifiable, the transition being completely smooth. 

We then solve the mean field equation numerically. All parameters are known, except 
for the width of the distribution of ball diameters. Diameter scatter is only qualified, in the 
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experiment, through a tolerance of ±4//m. We assume a uniform distribution whose width 

W is our single fitting parameter. We find that the best fit is obtained for W = 2.04 /xm, 

compatible with the tolerance figure. The fit itself is seen to be very good in the whole 

experimental range. The EMA is shown both for frictional and non-frictional balls. The 

results are only weakly sensitive to the value of rj, at striking variance with the sound 

velocity itself. For the latter quantity, however, the r\ dependence enters primarily through 

the G EFF /K EFF ratio. 

The experimental /-range is limited, on the low force side, to / ~ 2N. Furthermore, 

no direct measurement of the equation of state f(a) is available yet. It is therefore of 

great interest to complement the above comparison with a confrontation between mean field 

predictions and the simulations performed by Roux J7j on the same system for frictionless 

balls with a uniform random distribution of diameters. Applied forces span the whole range 

from almost zero up to the upper experimental limit. The comparison for the equation 

of state is shown on Figure [], which is drawn in the dimensionless representation used in 

reference J7J: 

d + lW-a „ f 



r 



W ' J \E*d*W* 
In these units, the exact equation of state f*(a*) is a universal function 0, if approximation 

( p2|) is used. It is easy to verify on Eq. fl38|) that the EMA has the same property. 

Clearly, a* > 0. For a* = 0, the balls barely start touching, a* = .343 is the so-called 

rigidity threshold, below which the disordered lattice cannot sustain compression, a* > 1 is 

the region where all neighbors are already in contact, and a* ^> 1 is the high pressure limit. 

In a narrow range above the rigidity threshold, a quasi-Hertz regime is found, followed, in 

the intermediate force range, by a steeper variation. The mean field result shown here, also 

calculated for frictionless balls, appears to agree with the numerical data for (a* — a*) ^ 0.3. 

Also shown are the curves corresponding to the ALA and to the AFA. They frame the exact 

and the EMA results from below and from above, respectively. For high pressures, all plots 

merge. We replot the same data in Fig. [7] using the linear scale for the lattice spacing, 
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and the | power scale for the forces, so that the natural range of various regimes can be 
better assessed. In particular, note the striking improvement of EMA over the AFA for 
intermediate pressures. 

In Fig. [8] we return to the experimentally more relevant stiffness-force dependence and 
plot in dimensionless form the mean-field and the ALA curves, the data from Roux's nu- 
merical simulations and their power law approximants, and also the experimental data 
of Gilles and Coste. These data are scaled using W = 2.04 /im and the force unit 
"^E*d^W^ = 19.89 N. The overall agreement is truly satisfactory, which is all the more 
remarkable that, while experiments and mean field are concerned with frictional balls, the 
simulations relate to the frictionless case. This is to be related with our previous observation 
that, in the h.c.p. lattice studied here, the normal stiffness is only weakly sensitive to shear 
interball forces. Such might not be the case for other ball arrays. 

As is well known, since mean field theories are not systematic expansions in powers of a 
small parameter, they do not allow for a precise direct assessment of their range of validity. 
From this discussion, we can state empirically that the validity of our effective medium 
approximation is limited to dimensionless forces /* ^ 0.08. This we confirm by calculating, 
within the mean field approximation itself, the relative force fluctuations Af*/f* = Af/f, 
and AF AVE /F AVE , where 

W) 2 = E c Q {FQ{a) - K^fl - /(a)) 2 ) (50) 

w 

The plot of Figure |9| shows that the relative fluctuations decrease rapidly with increas- 
ing displacement, i.e. pressure, the above mentioned empirical limit corresponding to the 
reasonable value Af*/f* ~ 0.6. The reduction of the local stress (force) fluctuations due to 
the self-consistent local displacements is reflected in the relative magnitude of the EMA and 
AFA fluctuations. The frictionless balls appear to adjust their positions more effectively, as 
could be expected. 

Finally, we have estimated the average fraction of active contacts per ball N c in our 
effective lattice in the following way. We define it to be the average number of neighbors 
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with an intercenter distance smaller than the sum of the corresponding radii. 

N c = £ c Q c Q , tf (^±^ - a -u Q ) (51) 

Q,Q' 

Note that this expression does not derive systematically from the mean field formalism, 

in which the notion of contact number does not enter explicitly. It should therefore be 

considered as indicative only. N c is plotted on Figure |l^ together with the values calculated 

by Roux. Expression ([y]) is seen to systematically underestimate N c . It appears that the 

validity of the mean field extends down to N c ~ 85%. However, it is important to point out 

that the pressure range where the sound velocity departs from the Hertz behavior extends 

well into the pressure range where the contact fraction has already saturated to 1. This 

corroborates strongly the idea that it is not the connectivity itself which is responsible for 

the non-Hertz behavior, but the presence of disorder-induced stress fluctuations. This is 



shown in Fig. [11], where we plot together the EMA results for N c , Af*/f* and the ratio 
k/k ALA . This latter quantity is a direct measure of the deviations of the effective stiffness 
from the ideal Hertz law. No sharp change of regime occurs at saturation of iV c , the non- 
Hertz behavior and stress fluctuations appear to extend to higher pressures and gradually 
tend to zero together. 

V. CONCLUSION 

On the basis of the above discussion, we are therefore able to conclude that our mean 
field theory provides quite a satisfactory description of the pressure dependence of the bulk 
mechanical properties, and hence of the sound velocity, in an array of frictional balls in the 
"high" pressure range. This corresponds to the regime in which the ball network is strongly 
overdetermined, that is where connectivity is close to its saturation value. 

We believe that this agreement, which permits to account for the non-Hertzian behavior 
down to the v « 1/4 range, is due to the fact that this theory does capture, though in an 
approximate manner, the existence of disorder induced stress fluctuations, and that these 
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are self-consistently related with the global mechanical state of the system. In other words, 
the true disorder strength in the problem is not intrinsically given by the dispersion of 
unstressed ball diameters, but determined from this together with the elastic deformation 
field. 

We have based our single-site description upon the simplest possible approximation, 
which amounts to a spherical averaging of local lattice distorsions. That this is sufficient to 
produce a satisfactory theory must certainly be attributed to the high connectivity of the 
h.c.p. structure. It should for this reason be applicable as well to a 3D lattice such as the 
fee one. 

Systematic improvement upon this approximation is possible. This should be based on 
extending the basic building block from our (central ball + average bond star) to a full cluster 
made of the central ball and of its six neighbors. Such an extension, although obviously 
heavy, would have the merit of explicitly allowing for local asymmetric configurations, which 
we have ignored here. It would also be a first step towards taking into account spatial stress 
correlations, which are completely overlooked in the present approach. 

In view of this last remark, the good agreement obtained here with sound velocity data 
calls for a physical comment. Coherent sound measures a mechanical response on the scale of 
the corresponding effective wavelength which, in the low frequency regime of the Gilles-Coste 
experiments, is at least of the order of 10 ball diameters, as discussed in Sec. IV. It is now 



well documented | 18H , [ fLTjj , [21 that the correlation length of the stress network in granular 
packings is, at most, of order a few d. It is therefore to be expected that mean field theories, 
though by nature unable to capture long range correlations, are well adapted to describe 
large scale properties. That is, when focussing on large scale mechanical responses, the 
pertinent notion is that of a global fluctuating stress network. The complementary notion of 
stress chains, which emphasizes the long range part of correlations, is certainly more relevant 
to the question of acoustic scattering, which becomes essential at higher frequencies. 

Obviously, an important pending question is concerned with the possibility of extending 
the mean field approach to the more important case of topologically disordered random 
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grain packings. Such a step, however desirable it might be, is by no means straightforward, 
as can be inferred from the, yet insuperable, difficulties which have been met when trying 
to extend the theories of electronic and vibrational properties of random substitutional 
alloys to amorphous materials. These can be assigned to the absence of a natural reference 
configuration possessing long range order. Up to now, the only existing theoretical frame for 
amorphous solids is a structureless, homogeneous effective medium. This would amount in 
the present problem to treating the effective medium in the continuum limit, which is clearly 
inedequate to account for stress fluctuation effects. In this perspective, numerical studies 
appear essential as a basis for trying to build up the needed original theoretical concepts. 

However, we believe that the qualitative idea which emerges from this work, namely that 
it is the disorder induced stress fluctuations which are responsible for the pressure depen- 
dence of the sound velocity in granular packings, will carry over, as well, to topologically 
disordered systems. 
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APPENDIX A: CALCULATION OF a, (3 

We want here to calculate the expression of the coefficients a, (3 defined in Section III. 
For this purpose, we need to solve the following problem. Consider an ideal h.c.p. lattice, 
with interball distance at equilibrium d, in which we singularize a central site (0) . Apply to 
the six balls i — 1, . . . , 6 of the corona of its nearest neighbors excess forces directed along 
the (Oi) bonds F; = Fhi (see Figure [|). Call Uj = u^i the resulting displacement of nearest 
neighbor (i). Then, by definition, a = Ui a /ui, where Uj Q is the displacement of the second 
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neighbor (ia) along the n^ direction; (3 = Uip/ui, with u^ the displacement of the next 
nearest neighbor along the direction of -73 n^ + n i+1 . 

In order to avoid excessive algebraic heaviness, and in view of the fact that numerical 
studies lead us to conclude to the very weak sensitivity of the mean field results to moderate 
variations of a, (3, we limit ourselves to the simple case of vanishing Mindlin shear stiffness: 
77 = 0. 

In this case rotations become irrelevant and, in the linear response regime, the force on 
a ball is related to the displacements of its neighbors by: 

f^y.%-^ ( A1 ) 

{3} 
So that, inverting in Fourier space : 

u l = ]Te^ R *5(q)F(q) (A2) 

k 

with 

5(q)=^ 1 (q) (A3) 

From equation(H), D(q) is a 2 x 2 matrix acting in (x, y) space, defined by: 



D(q) = E e iq -( R '- R j) Dij 

3 

= 2k J2 (1 — cos(dq.iip)) h p h p 

P =i 



(A4) 



The applied forces are defined by: 

F i = F [ni(5ii - S u ) + n 2 (^ 3 - 6 i6 ) + n 3 (fe - 5 i2 )] (A5) 

where the unit vectors h p and the sites (i) are labelled according to Figure |1|. 
Then: 

f 3 

F (q) = 4? E ("20 h p sin(rfq.n p ) (A6) 

iV p=i 

N being the number of sites in the lattice. 
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Taking advantage of the fact that, by symmetry, ni.(ui — 114) = 62.(113 — Uq) = 63.(115 
u 2 ), one gets: 



ni.^ui -u 4 J 
Ml= 2 



^EE (A7) 



(A8) 



^ q PfP=1 

n p .A(q).n r (cos[<iq.(fi p — n r )] — cos[<iq.(n p + n r )]) 
Analogously, the displacement U\ a of the second neighbor along Ox: 

^=^££ 

q p,r=l 

n p .A(q).n r (cos[<iq.(2n p — n r )] — cos[<iq.(2n p + n r )]) 

From equations (0),(fD, 

P(q) = detp(q)] = £ T p T p+1 (A9) 

p=i 

with 

r p = (1 - cos(rfq.n p ) r p+3 = r p (A10) 

Finally, one obtains, for a = u\ a ju\. 

E^(q) E r p (E p 2 i-Eg 2 )(EW 1 -E p 1 i) 
« = ^ ^ § (All) 

E^(q) E r p (££i - 4+ 2 ) 2 

q p=i 

where we have set: 

E p n) = sin(ndq.n p ) T p = sin(dq.t p ) (A12) 

A completely analogous calculation involving the displacement u±p of the second nearest 
neighbor yields: 
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/* = ^ ^-^ (A13) 

E^(q) E r^i-E^ 

q p=i 

Performing numerically the q-integrations over the Brillouin zone, we obtain: 

a = 0.585405 (3 = 0.232971 (A14) 
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FIGURES 

FIG. 1. Sketch of the ideal 2D h.c.p. lattice. 

FIG. 2. Dimensionless dispersion curves for vibrations in the h.c.p. lattice along two principal 
directions in the Brillouin zone. Units: ( k /M) 2 for frequency, 2n/d for wave vector. Upper panel: 
frictional balls, t/steel = 0.84. Lower panel: frictionless balls. Middle: Brillouin zone, TK = ^2n/d, 

TM = -^2-K/d. 

FIG. 3. Graphical representation of the EMA averaging process. On the left, a cluster (central 
ball surrounded by its corona, balls 1 -J- 6) embedded in the effective medium. Configuration 
averaging (...) restores the effective medium on the right. 

FIG. 4. Displacements of the neighbors of the cluster central ball 0. Corona balls: 1, 2, 6. 
Second neighbors: effective balls in the on-top (a) and bridge (f3 and /?') positions. 

FIG. 5. Square root of effective compressive stiffness vs. effective contact force. Squares: 
experiment (Ref. ||). Thick line: EMA for frictional balls. Dashed line: EMA for frictionless 
balls. Thin line: ALA (ideal Hertz dependence). 

FIG. 6. Effective force vs. displacement in dimensionless units (see text) for frictionless balls, 
a* corresponds to the rigidity threshold. Dots: numerical results. Dash-dotted lines: power law 
approximants of numerical results (Ref. ) . 

FIG. 7. Effective force vs. displacement in the "Hertz" representation. Data and symbols as 
Fig. g 
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FIG. 8. Square root of effective compressive stiffness vs. effective contact force (dimensionless 
units). 

FIG. 9. Relative fluctuation of local contact force vs. dimensionless displacement. Full line: 
EMA for frictional balls. Thick dashed line: EMA for frictionless balls. Dotted line: AFA. 
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FIG. 10. Active contact fraction vs. dimensionless displacement. Full line: EMA for frictional 
balls. Thick dashed line: EMA for frictionless balls. Dots: numerical simulation |7J. 

FIG. 11. The relative departure k* /k AhA of the EMA effective normal stiffness from ideal Hertz 
and the EMA force fluctuation Af* / f* exhibit no sharp change at the point /* ~ 0.3 where the 
active contact fraction N r saturates to 1. 



35 



6 
5 
4 
3 
2 
1 




" 'H steel 


'■■••N» 
'"■•• \ 


r^TSvX. 


.**>^ .■■''^""" 






V 


y^ 



M 



K 




5000 



- 1 1 r - 



- 1 1 r - 



4000 



C\J 



3000 



2000 



W = 2.04jLim 




10 



f[N] -> 



0.1 : 



0.01 7 



0.001 r 



0.0001 



1e-05 VjL 
0.001 




0.01 



0.1 



a - a 







CO 
* 




1 - 





1 ' ' 





1 1 1 1 ,.. 1 1 1 


B • 1/ i i i i 


VkVk 

I I 1 I I 




- 


• 




experiment 
simulation 

■ 1 



0.001 



0.01 



0.1 



10 



0.8 - 



0.6 



0.4 - 



0.2 












1 - 



£ 0.8 - 

O 

03 



O 
03 



O 
O 

CD 
> 

o 

03 



0.6 - 



0.4 - 



0.2 - 







- .._^^« , 

• / 



0.2 0.4 0.6 * 0.8 1 1.2 1.4 




- 1 1 1 1- 



0.6 



0.4 



0.2 







Af/f 



_i i i i_ 



0.1 



